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Abstract 

We consider the development of exponential methods for the robust time discretization 
of space inhomogeneous Boltzmann equations in stiff regimes. Compared to the space homo- 
geneous case, or more in general to the case of splitting based methods, studied in Dimarco 
Pareschi ;6] a major difficulty is that the local Maxwellian equilibrium state is not constant in 
a time step and thus needs a proper numerical treatment. We show how to derive asymptotic 
preserving (AP) schemes of arbitrary order and in particular using the Shu-Osher represen- 
tation of Runge-Kutta methods we explore the monotonicity properties of such schemes, like 
strong stability preserving (SSP) and positivity preserving. Several numerical results confirm 
our analysis. 

Keywords: Exponential Runge-Kutta methods, stiff equations, Boltzmann equation, fluid 
limits, asymptotic preserving schemes, strong stability preserving schemes. 
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1 Introduction 

The time discretization of kinetic equations in stiff regimes represents a computational challenge 
in the construction of numerical methods. In fact, in regimes close to the fluid-dynamic limit 
the collisional scale becomes dominant over the transport of particles and forces the numerical 
methods to operate with time discretization steps of the order of the Knudsen number. On the 
other hand the use of implicit integration techniques presents considerable limitations in most 
applications since the cost required for the inversion of the collisional operator is prohibitive 
therefore limiting such techniques to simple linear operators. 

In recent years there has been a remarkable development of numerical techniques specifically 
designed for such situations [H El [TOj, (H El [161 HO] • The basic idea common to these techniques is 
to avoid the resolution of small time scales by using some a priori knowledge on the asymptotic 
behavior of the kinetic equation. In particular, we recall among the different possible approaches 
domain decomposition strategies and hybrid methods at different levels [191 HI El El EZ] ■ 

Asymptotic-preserving schemes have been particularly successful in the construction of uncon- 
ditionally stable time discretization methods that avoids the inversion of the collision operator. 
For a nice survey on asymptotic-preserving scheme for various kinds of systems see, for example, 
the review paper by Shi Jin |15j . In the case of Boltzmann kinetic equations we also refer to the 
recent review by Pareschi and Russo |21j . 

In this paper we propose a new class of exponential integrators for the inhomogeneous Boltz- 
mann equation and related kinetic equations which is based on explicit exponential Runge-Kutta 
methods |14|, [T7] . More precisely we extend the method recently presented by one of the authors 
for homogeneous Boltzmann equations [6j to the inhomogeneous case by avoiding splitting tech- 
niques. The main feature of the approach here proposed is that it works uniformly with very 
high-order for a wide range of Knudsen numbers and avoids the solution of nonlinear systems 
of equations even in stiff regimes. Compared to penalized Implicit-Explicit (IMEX) techniques 
[H E] the main advantage of the class of methods here presented is the capability to easily achieve 
high order accuracy, asymptotic preservation and monotonicity of the numerical solution. 
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At variance with the approach presented in Dimarco, Pareschi [6] here we used the Shu-Osher 
representation of Runge-Kutta methods |26| . This turns out to be essential in order to obtain 
non splitting schemes with better monotonicity properties (usually referred to as strong stability 
properties [12]), which permits for example to obtain positivity preserving schemes. In particular 
we construct methods which are uniformly accurate using two different strategies. The first class 
of methods is based on the use of a suitable time independent equilibrium state which permits to 
recover high order accuracy and positivity of the numerical solution. However since the method is 
based on a constant equilibrium computed at the end time it may suffer of accuracy deterioration 
in intermediate regimes. The second class of methods is based on computing explicitly the time 
variation of the Maxwellian state. This permits to obtain schemes with better uniform accuracy 
but loosing some of the monotonicity property obtained with the first technique. 

The rest of the manuscript is organized as follows. In the next section we introduce some pre- 
liminary material concerning the Bolztmann equation and its fluid-limit. In Section 3 we derive 
the novel asymptotic-preserving exponential Runge-Kutta schemes. Two different approaches 
are presented. The properties of the two approaches are then studied in Section 4. In partic- 
ular monotonicity properties are investigated. Finally in Section 5 several numerical results for 
schemes up to third order are presented which show the uniform high order accuracy properties 
of the present methods. Some theoretical proofs are reported in a separate appendix. 



2 The Boltzmann Equation and its fluid-dynamic limit 

2.1 Boltzmann Equation 

The Boltzmann equation describes the evolution of the density distribution of rarefied gases. 
We use f(t,x,v) to represent the distribution function at time t on the phase space (x,v). The 
Boltzmann equation is given by 

dtf + v V x f = - e Q{f, /), t > 0, (x, v)€R d x R d , (2.1) 

with 

Q(f,f) = Q + -fQ~ = [ I (f'fl-ff*)B(\v-v*\,u)dv*(k>. (2.2) 

Here, B is the collision kernel, e > is the Knudsen number, uj is a unit vector, and S d ~^ is 
the unit sphere defined in R d space. We use the shorthands /' = f(t,x,v') and = 
There are many variations for the collision kernel B. One simple case is the case of Maxwell 
molecules when 

b = b( 9 ^ 

V \9\ 

with the relative velocity g = v — v* . 
The collisional velocities v' and v'^ satisfy 

v' = v- -(g - \g\u), (2.3a) 

vl = v* + -(g-\g\u). (2.3b) 
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This deduction is based on momentum and energy conservations 

v + v * = v' + v[ , 

I 1 2 i I |2 I / 1 2 i I / |2 

hi + Kl = \ v I + Kl • 

In ti-dimensional space, we define the following macroscopic quantities p is the mass density 
(here we assume mass is 1, thus number density and mass density have the same value); u is a 
(i-dimensional vector that represent the average velocity; E is the total energy; e is the specific 
internal energy; T is the temperature; S is the stress tensor; and q is the heat flux vector, given 
by 



fdv, pu = J vfdv, 

E=-pu 2 + pe = -J \v\ 2 fdv, e = -T = -Jf\v-u\ 2 dv J (2.4) 

S = J (v — u) ® (v — u)fdv, q = - J (v — u)\v — u\ 2 fdv. 

2.2 Conservations and fluid limit 

Cross section may vary, but the first d + 2 moments of the collision term are always zero. They 
are obtained by multiplying the collision term with (f> = (l,v, ^\v\ 2 ) T and then integrating with 
respect to v, i.e. 

<Q> = J Q(f)dv = 0, 
< V Q > = J vQ(f)dv = 0, 

<\v 2 Q> = j\\v\ 2 Q(f)dv = 0. (2.5) 

Based on these formulas, when taking moments of the Boltzmann equation, one obtains mass, 
momentum and energy conservation 

dtP + ^x ■ (pu) =< Q >= 0, 

d t (pu) + V X -(S + pu 2 ) = 1 < vQ >= 0, 

£ 

d t E + V X - (Eu + Su + q) = - < \\v\ 2 Q >= 0. 

For small values of e, the standard Chapman-Enskog expansion around the local Maxwellian 

/ J 1 Y /2 ( (v-u(t,x)) 2 \ , , 

M(t,x,v)= p(t, x) - — , exp - V — ^ , (2.6) 



2irT(t,x)J ^\ 2T(t,x) 
shows that at the leading order the moment system yields its Euler limit 

dtp + V • (pu) = 0, 

d t (pu) + V- (pu®u + pTI) = 0, (2.7) 
8 t E + V • ((E + pT)u) = 0, 
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where I is the identity matrix. 



3 Exponential Runge-Kutta (ExpRK) methods 



In this section we would like to extend the Exponential RK method in [6] for the homogeneous 



Boltzmann equation to the inhomogeneous case (2.1). It has been known for long that time 
splitting methods degenerate to first order accuracy in the fluid-limit (see [6] and the references 
therein) so, to achieve high order of accuracy in stiff regimes, time splitting should be avoided. 



3.1 Reformulation of the problem and notations 

To achieve AP property and robustness in stiff regimes, an implicit method should be adopted. 
However, due to the complexity and nonlocal property of the collision term Q, directly inverting 
it is prohibitively expensive. The Exponential Runge-Kutta method overcomes this difficulty by 
transforming the equation into the exponential form, and forces the solution to approach to the 
equilibrium that captures its asymptotic Euler limit as e tends to zero, thus it is an AP scheme. 
Following the approach in [6], one can define 



P = Q + fif, /U > 0. 



(3.1) 



Let us now consider a nonnegative function M, hereafter called the equilibrium function, and 



using (2.1) compute 



(f ~ M)e 



put Is 



d t (f - M)e^ t/£ + (/ - M)^e M ' A 



-(Q + lif-nM)-dtM-vV x f 



(3.2) 



(P _ M M) -d t M-v V x f 



Note that the equation above is equivalent to the original Boltzmann equation (2.1) as long as 
fi is independent on time. In the simplified case of the BGK collision operator Q = fi(M — /), 
where M is the local Maxwellian given by (2.6), the problem reformulation just described applies 
with P = [iM. Moreover there is no requirement on the form of M at all - it can be an arbitrary 
function. However, to obtain AP property, one has to be careful in picking up its definition, so 
that the correct asymptotic limit could be captured. 

We analyze two different approaches in the following two subsections, and adopt a suitable ex- 
plicit Runge-Kutta scheme to solve them. For readers' convenience, we firstly give the expression 
of the Runge-Kutta method used here. Given a large set of ODEs 



dty = F(t,y), 



(3.3) 



obtained for example using the method of lines from a given PDE, if data y n at time step t n is 
known, to compute for the value y n+1 at t n+1 = t n + h, a classical u-step explicit Runga-Kutta 
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scheme for equation (3.3) writes 



i-l 



Step i: 



J nM =y n + h J2a ij F(t n + c j h,y n ' 



Final step: y n+l = y n + h ^ biF{t n + ah, y n ' {i) ) , 



(3.4) 



where X)}=i a «i = c *' Si ^ = 1> an d stands for the estimate of y at i = t n + q/i. Different 
Runge-Kutta method gives different set of coefficients. In the sequel we drop superscript n for 
evaluation of y at sub-stages and use = y ra >W . 

Another form of RK method which has proved to be useful in the analysis of the mono- 
tonicity properties of Runge-Kutta schemes is the so-called Shu-Osher representation [26J. This 
representation is essential in the study of the positivity properties that will be carried out later 



i-l 



Step i: 
Final step: 



,n+l 



^[^y U) + h^F(t n + Cj h,y^) 



3=1 
v 



£ [a u+lj y^ + hf3 u+lj F(t n + Cj h,y^) 



(3.5) 



j=i 



Let us point out that this latter representation is not unique. Here a.%j are parameters such that 
S}=i a ij = Without loss of generality, it is natural to set 



Pij — OC%j (Cj Cj) 



3 1 ' 



(3.6) 



for consistency. 



Remark 1. Expression (3.6) is equivalent with the classical one which says [26] 

i-l 

Pij = a ij ~ ^2 aikCLk 3- ( 3 - 7 ) 
k=j+l 

In fact, assume one has yW = y n + h Y^k=i a jkF <yk \ V j < i, where F^ is a shorthand for 
F(t n + c k h,yW), then, by (3.6) one has 



,0 



i-l 

n 

3=1 

£ 

j<i 



otijy 



Cj) 



aij(ci - Cj)hFV) 



a 



l 3 



y n + h^2a ]k F^ ] + onjia - Cj)hF® 
k<j 



i-l 



y n + h ^ a^a/cj + ctijici - cj) 

3<i \k=j+l 



F (3) 



(3. 



This clearly requires a%j = £%(cj — cj) + ^ a^a^j. Given (3.6), it is aij = fyj + ^ ctik a kj> which 
is exactly the classical Shu-Osher representation. 
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3.2 Exponential RK schemes with fixed equilibrium function 



Since the choice of the equilibrium function M in (3.2) is arbitrary, in this subsection, we assume 



M as a function independent of time in each time step, i.e. M is a function given a-priori. Thus 



(3.2) could be rewritten as 

d t [(/ - M)e^ t/£ 



: {P-fj,M)-v-V x f 



(3.9) 



Analytically, the equation (3.9) is equivalent to the original inhomogeneous Boltzmann equation 



as long as M is a function independent of time and [i is a constant. But the associated numerical 
scheme can preserve asymptotic limit only if M is chosen in a correct way, as will be clearer later. 
On the other hand fi plays a role in order to guarantee positivity of the numerical solution as 
will be seen in section. 

Remark 2. Obviously M is required not to change in each time step. But for different time 
steps, we are free to use different functions. This is in fact what we will do, we evolve M before 
each time step with a suitable scheme, and then use this computed value function to construct the 
AP exponential scheme. 

3.2.1 The numerical scheme: ExpRK-F 



on the right of 
Step i: 



Compared to (SJi^y turns out to be (/ — M)e Mi//£ and the associated evolution function F(t,y) 

-(P — jtiM) — v ■ V x f e^*/ £ . Thus we have the following scheme 



3.3) is 



i-l 



(/W - M)e c * A = (f n - M) + <Hj~ P {j) -fiM-ev- V x f ij) 



Final Step: (/ n+1 - M)e A = (/**- Af) + &»- \p {i) - /uAf - ev ■ V x f 



i=i 



r(i) 



e°i x , 



Ci\ 



where we used A = and = P(f^) for simplicity. Simple algebra gives 
• Step i: 



(3.10) 



i-l 



i-l 



f 



(i) 



l _ e -ciX _ W..v As --' 

3=1 



^o ii Ae A( - c ' +c ^ ] M+e- c ^f n +^a tj \e^ 

3=1 



A(c i -c i ) 



'i 



Final Step: 



/ n+1 = ^1 - e" A - ftiAe A (- 1+Ci M M + e- x f n + ^ Me A(ci_1) ( 
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3.2.2 Choice and evaluation of M 



If it is assumed that 



= C\ < C2 < . . . < c„ < 1, 



(3.11) 



then the same arguments used in [B] shows immediately, that as e — > 0, A — > oo, the scheme 
pushes f n+1 going to M. So to obtain AP property, M above should be the Maxwellian at time 
level n + 1 that has the right moments. To get the right moments, the simplest way is to evolve 
the corresponding macroscopic limit equations, say the Euler equation. We propose solving the 
Euler equation first to obtain the macroscopic quantities of the Maxwellian for the next time step, 
and make use of them to define M. To achieve high order for all regimes, both the macro-solver 
and micro-solver should be handled by numerical schemes with the same order of accuracy in 
space and time. The most natural way in time discretization is the explicit Runge-Kutta scheme 
using the same coefficients as the one for the kinetic equation 



Step i : 



Final Step: 




n+l 




( pu 

pu® u + pT 
V (E + P T)u 
I pu 

pu®u + pT 
\ (E + P T)u 



(?) 



(3.12) 



Remark 3. 



Note that this method gives us a simple way to couple macro-solver with micro-solver. When 
e is considerably big, the accuracy of the method is controlled by the micro-solver. And as 
e vanishes, the method pushes f going to M, which is defined by macroscopic quantities 
computed through the Euler equation while the order of accuracy is given by the macro- 
solver. 

In principle it is possible to adopt other strategies to compute a more accurate time indepen- 
dent equilibrium function in intermediate regions. For example one can use the ES — BGK 
Maxwellian |2J/ at time n + 1 or one can use the Navier-Stokes equation as the macro- 
counterpart. Here however we do not explore further in these directions. 



The assumption (3.11), although strongly simplifies computations, is in fact not necessary 
to prove asymptotic preservation. In fact such assumption is independent of the structure 
of the operator P(f, /). We refer to Section 4-2 for more details. 



3.3 Exponential Runge-Kutta schemes with time varying equilibrium function 

The approach just described has the nice feature of being extremely simple to construct and 
implement. As we will see in the next section it also possesses several nice features concerning 
monotonicity. On the other hand it is clear that choosing the limiting equilibrium state in the 
construction may produce a lack of accuracy in intermediate regimes. To overcome this aspect 
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here we consider the most natural choice of equilibrium function, namely the local Maxwellian 
equilibrium state M = M. The major difficulty in this case is due to the time dependent nature 
of such equilibrium function. 
Now rewrite the equation as 



(/ - M) exp 



/.d 



P-fiM 



v-V x f- d t M exp 



(3.13) 



and here we define M has a Gaussian profile that shares the same first d + 2 moments with /. 
The moments' equations are governed by 



d t / 4>fdv 



bv ■ V x fdv = 0, 



(3.14) 



with 



1 v — 

5 U 1 2 



3.3.1 The numerical scheme: ExpRK-V 

The Runge-Kutta method is adopted to solve the system 

' d t {f - M)e^ t/e = 1 (P -\xM-ev- V x f - £d t M)e> lt/£ , 



d t / 4>fdv 



J 4>v ■ V x fdv. 



Thus we have the following scheme 
Step i: 



i-1 



(/« _ M«)e c * A = (f n - M n ) + aij - - fihl^ - ev ■ V x f {j) - ed t M^ 



f(/>f®di 



Final Step: 

(p+l - M n+1 )e x 
f (j)f n+1 dv 



J ^n dv + ^ aij f-h I cf>v ■ V x f^dv) ; 

.7=1 ' 



(3.15a) 



(f n - M n ) + h- P [i) ~ - ev ■ V x f {i) - ed t M® 

i=i £ 

J 4>f n dv + Yj>i [~h J <f>v ■ V x f®dv) . 



(3.15b) 



The first equation in (3.15a) shows that in each sub-stage i, to compute for besides the 
known /w and easily obtained M®, one also needs dtM^\ P®, v ■ V x /^ for all j < i, and 
that is evaluated at the current time sub-stage. 
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3.3.2 Computation of M and d t M 

Here we show how to compute M« and dtM^\ 

Computation of M^> : 



solve the second equation of (3.15a), to get evaluation of macroscopic quantities at t n + Cih. 
Then the Maxwellian M' 1 ' is given by (2.6). 

Computation of dtM^ : 
Write d t M as 

d t M = d p Md t p + V U M • d t u + d T Md t T, (3.16) 
and dtp, dtu and dtT can be computed from taking moments of the original equation 



t 



pu 



Ot 



v | Mdv = d t 

2 
1 

v I v ■ V x fdv. 



1 

v I fdv 

2 



(3.17) 



To be specific, with data at sub-stage (J) in d-dimensional space, one has 
dtM® = d p M^d t p U) + V M M^ ■ d t u {j) + d T M^d t T^\ 

with 



dpM® 



,0') 



drM® = 



2(r0'))2 2rw 



and 



v ® v ■ V x f ij) dv , 



dpV) 



2E® 



d t p U) - 2p^u^d t u ij) - v 2 v V x f U) d 



(3.18) 
(3.19a) 

(3.19b) 
(3.19c) 

(3.19d) 



The dtp and dtu term in (3.19d) is evaluated by ( 3.19b[ ) and (3.19c). Clearly, all other 
macroscopic quantities p^\ vS^ and are associated to 
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4 Properties of ExpRK schemes 
4.1 Positivity and monotonicity properties 

Usually positivity, although very important for kinetic equations, is extremely hard to be obtained 



when using high order schemes. Here we show that thanks to the Shu-Osher representation (3.5) 
we can follow to prove positivity (and hence SSP property) for the fixed M method ExpRK-F. 
Before proving the theorem we make the following assumption. 

Assumption 1. For a given f > there exists h* > such that 

f -hv-V x f >0, V0</i</i*. 

The above assumption is the minimal requirement on / in order to obtain a non negative 
scheme. Next we can state 



Theorem 1. Let us consider an ExpRK-F method defined by (3.10), and fiij > in (3.6). Then 



there exist h* > and /x* > such that f n+1 > provided that f n > 0, fi > [i* and < h < h*. 
Proof. Using the Shu-Osher representation, one could rewrite the scheme as 

Step i: (/(*) - M)e c ' A = ^ e c > x jc*y(/ (i) — M) + ^ P U) - (J,M — ev ■ V x f ( - j) j 

Final Step: (/ n+1 - M)e x = ^ e c J x ia v+Xj {P — M) + Pu+ij^ P (j) - f*M - ev ■ V x f^ j 



j 

Simple algebra gives, for V, i = 1, • ■ • ,v, j < i 



/« =M 1 - ^~ Ci)X + A/%) 



i-l 



(cj- Ci )A 



pO') 



i=i 

i-i 



+ £ a,,e^- c ^ f /W - ^ • V„/^ . (4.1) 

The same derivation can be also carried out for the final step. If this is a convex combination, 
then, to have positivity, one only check that each of them is positive 

M > 0; (4.2a) 

P {j) > 0; (4.2b) 

f U) _ h ^l v . Vxf U) > o. (4.2c) 

Positivity of M is obvious, and P^' is positive if one has big enough /j, 

\i > fi* = sup \Q~\ => P = Q + /// = Q + - fQ~ +nf>0. 
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To handle (4.2c), one just need to adopt Assumption 1. It is positive if 

• ( °ij 
mn I — jL , 

ij \Pij 



< h < K = min -^-h* 



which guarantees (4.2c) 



To check the convexity of ( |4.1[ ), it should be proved that 

Y,^ iCj ~ CH)X (^ + XP l j)<l- (4.3) 

j 

This can be seen by just taking the derivative with respect to A. Use Aij = a — Cj 



d_ 



X/ A " A (n„ • A/%) 
i t 

e~ A ^ x (-Aijiaij + A/%) + /%■) (4.4) 

j 

e~ A ^ x {-pijAijX + /% - A^) < (4.5) 



In the last step, /%■ = ctij(ci — Cj) is used. Thus the left-hand side of expression (4.3) is mono- 
tonically decreasing with respect to A > and has a maximum 

E a v = 1 



at A = 0. Similarly we can proceed for the final step. This confirms (4.3) and finishes our 



proof. □ 

Since the proof above is based on a convexity argument, we also have monotonicity of the 
numerical solution or SSP property. Thus the building block of our exponential schemes is 
naturally given by the optimal SSP schemes which minimize the stability restriction on the time 
stepping. We refer to |12| for a review on SSP Runge-Kutta schemes. 

Remark 4. 

• Note that the proof above does not rely on the value A take, i.e. the scheme is positive 
uniformly in e. For the choice of fi* we refer the reader to the discussion in J3|> liOj/ . 

• Optimal second and third order SSP explicit Runge-Kutta methods such that (3ij > have 
been developed in the literature. However the classical third order SSP method by Shu and 
Osher \2b^ does not satisfy Cj < q for j < i. Note that standard second order midpoint and 
third order Heun methods satisfy the assumptions of Theorem 1 (see Table 1.1 page 135 in 

m)- 

• In it was proved that allfour stage, fourth order RK methods with positive CFL coeffi- 
cient h* must have at least one negative fy. The most popular fourth order method using 
five stage with nonnegative 0ij has been developed in 1241 ■ In [2^ the authors also proved 
that any method of order greater then four will have negative fiij. 
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• Positivity of ExpRK- V schemes is much more difficult to achieve because of the involvement 
of the dfM term. However, we can prove: 

(i) p is positive; 

(ii) the negative part of T is 0{he). 

We leave both the proofs of the above results to the appendix. 



4.2 Contraction and Asymptotic Preservation 

In this section, it will be presented that the new exponential Runge-Kutta schemes preserve 
the asymptotic limit of the Boltzmann equation. The proof is done by following the proof of 
contraction in [6]. 

If one check the formula (3.10) and (3.15b), it seems clear that under assumptions (3.11) 
the big A on the shoulder of exponential will push the distance between / and the Maxwellian 
function going to zero. But sometimes the Runge-Kutta method may have tough coefficients, say 
c v = 1. When this happens, the argument cannot be carried through. However, one could still 
prove AP property using the particular structure of the collision operator following the framework 
below. 

We need to make use of the following assumption. 

Assumption 2. There is a constant C big enough, such that \P(f,f) — P(g,g)\ < C\f — g\ 
where \-\ denotes a proper metric. 

Part of the proof for the metric d2 defined in P s (JR. d ) space (see [28] ) can be found in the 
appendix. 

Under this assumption, considering P(M, M) = Q(M, M) + \xM = jiM, one has 



\P(f,f)-liM\ <C\f-M\ 



(4.6) 



The derivation and the proof for both approaches being AP will be presented below. We first 
show that ExpRK-F is AP for any given explicit Runge-Kutta scheme. 

For AP property, one needs to show that as e — > 0, the scheme gives correct Euler limit. To 
do this, basically one needs to prove that / goes to the Maxwellian function whose macroscopic 
quantities solve the Euler equation (2.7). 

Let us define 



di 



v ■ V x f 



(0 



d 



f n -M 



/ w - M 

d= [d u d 2 ,--- ,d u ], D = [D 1 ,D 2 ,- - ,A, 

Moreover A is a lower-triangular matrix and E is a diagonal matrix given by 

A 



1,1, 



(4.7) 



A 



aije 



(cj-Ci)X 



E = diagje 



-ciA -C2A 



-e„A 



}• 



Lemma 1. Based on the definitions above, for ExpRK-F one has 

d< d (I - CA)' 1 . E • e + e (I - CA) _1 • A • D 
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Proof. It is just direct derivation from (3.10) 



i-l 



A 



(/« - M)e c * A = (/ n - M) + aij-e c ^(P ij) - fiM - ev ■ V x f ij) ) (4.8) 
3y taking the norm, adopting the triangle inequality, and make use of the assumption that 



P(f) - fiM 



< C 



f-M 



one gets 



A 



f(i) — m < f n -M e~ CiX + ^ Oij-^ Ci ~ Ci)x [C f {j) - M 



+ e 



• VJ« ) (4.9) 



Written in the matrix form, it becomes 

/ di \ ( d \ 

d 



d-2 
\d u J 



< E 



+ CA 



\ do j 



( di \ 

d 2 

\d u J 



+ eA 



D 2 
\ Du J 



Thus 

d< d E ■ e + CA ■ d + eA ■ D 

d< d (I - CA)' 1 • E • e + e (I — CA) -1 • A • D 

which completes the proof. 
Lemma 2. Define 

R^X) = e ~ x (l + ^6 • E" 1 (I - CA)' 1 E • Sj 

R 2 (\) = — e - x b ■ E^ 1 • (I - CAy 1 ■ (I + CA) 
ft 



then for scheme (3.10) we have 



< 



f n -M 



R 1 (X) + R 2 -D 



Proof. It is just a simple derivation. Define 

h 



k, = - {P^ - fiM - ev ■ V x f {i) )e 



Evidently, the previous lemma leads to 



|A:| < — E ■ (Cd + eD 
ft 



(4.10) 
(4.11) 

□ 



(4.12) 
(4.13) 

(4.14) 

(4.15) 

(4.16) 
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Back to (3.10), one has 



which implies 

fn+l _ M 



en+1 



X 



M 



f n -M) e- x + ^biki 



3=1 



<d e~ x + -e~ x P • E" 1 • (Cd + eD 
H V 



<e~ x [do + ^b-E- 1 ■ [(' (i- CAT' ■ ( d (] E ■ <■'+ :A ■ D ) + :D 



A: 



<d$e 



-x 



CA; 



1 + — b ■ E • (I - CA)" 1 • E • e 



sX 



+ —e~ x b ■ E" 1 • (I - CAY 1 ■ (I + CA) ■ D. 



(4.17) 

(4.18a) 
(4.18b) 
(4.18c) 
(4.18d) 



Here b = [bi, 62, • • • , b v ] is a row vector. The result (4.11) is also used. Plug in the definition of 
Ri and R2, one gets 



n+l 



M 



< 



f 



M 



i?i(A) +R 2 (X) ■ D. 



(4.19) 
□ 



The two lemmas above gives us the estimation of the convergence rate towards the Maxwellian. 
The smaller R\ is, the faster the function converges. R2 represents the drift from the transporta- 
tion, and is expected to be small in the limit. Also, the matrix A is usually a lower triangular 
matrix, and a strict lower triangular matrix for explicit Runge-Kutta, thus it is a nilpotent. 



Theorem 2. The method ExpRK-F defined by (3.10) is AP for general explicit Runge-Kutta 
method with < c\ < C2 < • • • < c v < 1 . 

Proof. Obviously if R\(X) = 0{e) and i?2(A) = 0(e) for e small enough, the theorem holds. In 
fact, for explicit Runge-Kutta method, A is a strict lower triangular matrix, and thus a nilpotent, 
then one has the following 



E" 1 (I - CA) -1 E = E" 1 (I + CA + C 2 A 2 + ■■■ + C l/ " 1 A I/ - 1 ) E 



I + B + B 2 H hB^ 1 



(4.20a) 
(4.20b) 



where A u = 0, definition B = CE -1 AE and E- X A 2 E = E^AEE^AE are used. According to the 
definition of A and E, it can be computed that 

CA 



CA„ ( ' < A ' A 



"13 ■ 



Thus I + X]fc® fc i s a ma trix such that: the element on the fcth diagonal is of order 0(X k ). This 



leads to obvious result 



Ri(\) 



l + ^-E-^l-Ar 1 !-^ = 0(e- x X y - 1 ) < 0(e) 
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Similar analysis can be carried to i?2(A) to show that it vanishes to zero as e — > 0. 
So as e —7- 0, |/ n+1 — M\ — > 0. By definition, M is defined by macroscopic quantities computed 
directly from the limit Euler equation, thus the numerical scheme is AP, which finishes the 
proof. □ 

The derivation of the scheme ExpRK-V is essentially the same, and in the end, one still has, 
in a condense form 

lyn+1 _ M n+1\ <\fn_ M n| R ^ + R 2 . £) ( 4 . 2 1) 



with Ri, R2, E, A defined in the same way as in kj ), but Di = \v V x f® + dtM®\. Following 
the same computations, one could prove that this method is AP too, but the proof is omitted for 
brevity. 



Theorem 3. The method ExpRK-V defined by (3.15) is AP for general explicit Runge-Kutta 
method. 



5 Numerical Example 

5.1 Convergence Rate Test 

In this example, we use smooth data to check the convergence rate of both methods. The problem 
is adopted from [8|: 1 dimensional in x and 2 dimensional in v. Initial distribution is given by 

f(t = 0,x,v) = ^j 1 f e T oW + e T oW J (5.1) 

with 

p (x) = ^(2 + sin(2 TO )), 

ui(x) = [0.75,-0.75] T , u 2 (x) = [-0.75, 0.75] T , 
T (x) = ^(5 + 2cos(27rx)). 

Domain is chosen as x G [0, 1] and periodic boundary condition on x is used. Note that the defi- 
nition of po, u\/2 and To do not represent the number density, average velocity and temperature. 

As one can see, the initial data is summation of two Gaussian functions centered at u\ and 
U2 respectively, and is far away from the Maxwellian. To check the convergence rate, we use 
N x = 128, 256, 512, 1024 grid points on x space, and N v = 32 points on v space. Time stepping 
At is chosen to satisfy CFL condition with CFL number being 0.5. We measure the L\ error of 
p and compute the decay rate through the following formula [29J 

\\pAx(t) - P2Ax(t)\\l 
error a -r = max n — , (5.2) 

*=*" ||P2A*(t)lll 

with Ax = Theoretically, a fcth order numerical scheme should give error < C (Ax) k for 
Ax small enough. 
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We compute this problem using spectral method |18j in v, WENO of order 3/5 |25j for x. 
For time discretization, we use the second and third order Runge-Kutta from [13J, Table 1 page 
135. We denote the four schemes under consideration as ExpRK2-F, ExpRK2-V, ExpRK3-F and 
ExpRK3-V. 

We compute the problem using the Maxwellian, and a distribution function away from the 
Maxwellian given above as initial data, for e = 1,0.1, 10 ,10~ 6 . Results are shown in Figure 
5.1 We also give the convergence rate Table 5.1 One can see that in kinetic regime, when e = 1, 
the two methods are almost the same, but as e becomes smaller, in the intermediate regime, for 
example e = 0.1 for the second order schemes and e = 10~ 3 for second and third order schemes 
with Maxwellian data, ExpRK-V performs better then ExpRK-F. In the hydrodynamic regime, 
however, the two methods give similar results again shown by the two pictures for e = 10~ 6 . It 
is remarkable that the third order methods achieve almost order 5 (the maximum achievable by 
the WENO solver) in many regimes. 



Initial Distribution 


Maxwellian Initial 


Non-Maxwellian Initial 


N x 


128 - 256 - 512 


256 - 512 - 1024 


128 - 256 - 512 


256 - 512 - 1024 


e = 1 


ExpRK2-F 
ExpRK2-V 
ExpRK3-F 
ExpRK3-V 


1.91327 
2.41608 
4.99725 
5.02508 


1.99502 
2.02347 
4.35014 
4.40379 


1.84968 
2.67733 
5.12959 
5.13515 


1.98504 
2.05436 
4.76788 
4.79080 


e = 0.1 


ExpRK2-F 
ExpRK2-V 
ExpRK3-F 
ExpRK3-V 


1.98218 
2.41411 
5.07621 
5.02220 


1.99539 
2.02293 
2.94707 
4.39651 


1.97725 
2.56620 
5.49587 
5.13859 


1.99454 
2.05830 
3.00335 
4.79264 


e = 10~ 3 


ExpRK2-F 
ExpRK2-V 
ExpRK3-F 
ExpRK3-V 


1.23711 
2.02344 
2.36140 
3.86882 


1.64976 
1.85924 
2.69178 
3.03223 


1.43331 
1.47466 
2.55225 
2.59114 


1.73501 
1.75496 
2.78275 
2.80353 


e = 10~ 6 


ExpRK2-F 
ExpRK2-V 
ExpRK3-F 
ExpRK3-V 


2.56137 
2.56137 
5.08829 
5.08830 


2.04519 
2.04519 
4.56695 
4.56704 


2.56137 
2.56383 
5.08830 
4.91909 


2.04519 
2.04859 
4.56699 
3.80638 



Table 1: Convergence rate for ExpRK methods with different initial data, in different regimes. 



5.2 A Sod Problem 



This simple example is adopted from [29J to check accuracy and AP of the numerical methods. 
It is a Riemann problem, and the solution to the associated Euler limit is a Sod problem. 

J (p,u x ,u y ,T) = (1,0,0, 1), if x < 0; 

1 (p, u x , u y , T) = (1/8, 0, 0, 1/4), if x > 0; 



In Figure 5.2 (left), we show that when e = 0.01 is comparably big, both the two new method 
proposed here match with the numerical results given by explicit scheme with dense mesh. Here 
the reference is given by Forward Euler with Ax = 1/500 and h = 0.0001. In Figure 5.2 (right), 
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= 1, Non-Ma: 



ExpRK2-F 
ExpRK2-V 
ExpRK3-F 
ExpRK3-V 



ExpRK2-F 

— O— ExpRK2-V 

— *— ExpRK3-F 

A ExpRK3-V 



1.1, Non-Ma: 









__ _o 


__ & ■ -~~ ' 


J* 


a" 


- ■ - ExpRKE-F 
-O- ExpRK2-V 
— *— ExpRK3-F 
-A- ExpRK3-V 



Q£. -9.5 - 

i .„_ 



ExpRK2-F 
ExpRK2-V 
ExpRK3-F 
ExpRK3-V 



e=10" J , Non-Ma 



ExpRKE-F 
ExpRKE— ' 
ExpRK3-F 
ExpRK3-V 
—2.1 



■ — . — ■ ExpRK2-F 

-O- ExpRK2-V 

— *— ExpRK3-F 

—A— ExpRK3-V 



epsilon = 1 D , Ma: 



epsilon = 10 , Non-Ma: 





















-9 










-9 








-9.5 










-9.5 


















5^ -1 D 








S -10.5 










jjf -10.5 








-11 










-1 1 
















- ■ - ExpRK2-F 








■ — . — ■ ExpRK2-F 


-1 1.5 








-O- ExpRKZ-V 
— *— ExpRK3-F 


-1 1 .5 






-O- ExpRK2-V 
-m— ExpRK3-F 










—A— ExpRK3-V 








—A— ExpRK3-V 



Figure 5.1: Convergence rate test. In each picture, 4 lines are plotted: the lines with dots, 
circles, stars and triangles on them are given by results of ExpRK2-F, ExpRK2-V, ExpRK3-F 
and ExpRK3-V respectively. The left column is for Maxwellian initial data, and the right column 
is for initial data away from Maxwellian (5.1). Each row, from the top to the bottom, shows 
results of e = 1/0.1/1CT 3 /10~ 6 respectively. 
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AP property is shown: it is clear that for e = 10 6 , numerical results capture the Euler limit 
the Euler limit is computed by kinetic scheme |22j . All plots are given at time t = 0.2. 



5.3 Mixing Regime 



In this example [29], we show numerical results to a problem with mixing regime. This problem 
is difficult because e vary with respect to space. As what we do in the first example, we take 
identical data along one space direction, so it is ID in space but 2D in velocity. An accurate AP 
scheme should be able to handle all e with considerably coarse mesh. Domain is chosen to be 
x G [—0.5,0.5], with e defined by 



'e Q + 0.5 (tanh (6 - 20x) + tanh (6 + 20a;))) x < 0.2; 
e x > 0.2 



(5.3) 



where eo is 10 3 . So e raise up from 10 3 to 0(1), and suddenly drop back to 10 3 as shown in 



Figure (5.3). Initial data is the give as 



f(t = 0, 



x, v 



Po(x) 
4irT (x) 



\v-u Q (x)\ 



I"+"oO)I z 



e 2T (x) _|_ g 2T (x) 



(5.4) 



with 



Po(x) 

Uq(x) 
Tq(x) 



2+sin (2ttx+tt) 



COS (2-7TX + 7r) 





(5.5) 



3+cos (2-irx+n) 
4 



Periodic boundary condition on x is applied. 

We compute the problem using both methods proposed in this paper together with standard 
explicit Runge-Kutta 2 and 3 in time used as the underline methods in the exponential schemes. 



Results are plotted in Figure 5.4 The reference solution is computed with a very fine mesh in 
time. Both methods give excellent results simply taking a CFL condition of 0.5 whereas explicit 
methods are forced to operate on a time scale 1000 times smaller. In particular, ExpRK3-V 
performs well uniformly on e by giving a more accurate description of the shock profiles. 



6 Conclusions and future developments 

In this paper we have presented a general way to construct high-order time discretization methods 
for the Boltzmann equations in stiff regimes which avoid the inversion of the collision operator. 
The main advantages compared to other methods presented in the literature is the capability to 
achieve high order uniformly with respect to the small Knudsen number and to originate monotone 
schemes thanks to the exponential structure of the coefficients. The approach presented here can 
be extended in principle to several other integro-differential kinetic equations where it is possible 
to identify a linear operator which preserves the asymptotic behavior of the system. For example 
in the case of the Landau equation this would involve the computation of the exact flow of the 
linear part, i.e. a matrix exponential, in the construction of the schemes. We leave this possibility 
to future research. 
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Figure 5.2: Consistency and AP. Left column: e = 0.01. The solid line is given by explicit scheme 
with dense mesh, while dots and circles are given by ExpRK2-F and ExpRK2-V respectively, 
both with N x = 100. h = Ax/20 satisfies the CFL condition with CFL number being 0.5. Right 
column: For e = 10~ 6 , both methods capture the Euler limit. The solid line is given by the 
kinetic scheme for the Euler equation, while the dots and circles are given by ExpRK2-F and 
ExpRK2-V. They perform well in rarefaction, contact line and shock. 
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L-(X] 



1 - Kinetic Regime- 




Euler Regime 



-0.5 -0.4 -0.3 -0.2 -0.1 0.1 0.2 0.3 0.4 0.5 



Figure 5.3: Mixing Regime: e{x) 
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7 Appendix 

7.1 Positivity of the mass density in ExpRK-V 



Theorem 4. The method ExpRK-V defined by (3.15) gives positive p, and the negative part of 
T is at most of order 0(he). 

To prove this theorem, we firstly check the following lemma. 

Lemma 3. In each sub-stage, the distribution function /w and have the same first d + 2 
moments. 

Proof. We prove this for sub-stage i. Assume for \/j<i, one has 



v (f^-M {j) )dv = 0. 



(7.1) 



V— 
2 



21 



t=0.75 



t=0.75 




Figure 5.4: The left column shows comparison of RK2 and RK3 using the ExpRK-V. The solid 
line is the reference solution with a very fine mesh in time and Ax = 0.005, the dash line is 
given by RK3 and the dotted line is given by RK2, both with N x = 50 points. The right column 
compare two methods, both given by RK3, with the reference. The dash line is given by ExpRK- 
V, and the dotted line is given by ExpRK-F. N x = 50 for both, h is chosen to satisfy CFL 
condition, in our case, the CFL number is chosen to be 0.5. 



22 



Then, one could take moments of the first equation in the scheme (3.15a), and gets 
(/W - M^)e c * x dv -- 



v 

2 



(f n - M n ) dv 




P^ - pM®) dv 



(7.2a) 
(7.2b) 



(ev ■ V x f ij) - ed t M( j A dv (7.2c) 



(7.2a) is zero for s ure, ( ] 7.2b ) is zero by definition of P and (7.1). (7.2c) is zero because of the 
computation from (3.18). Thus it is obvious that /W and share the same moments on each 
stage. □ 

With the previous lemma in hand, one could prove Theorem 4. 

Proof. As in the previous lemma, we only do the proof for sub-stage i. The final step can be 
dealt with in the same way. Rewrite the second equation of (3.15a) in Shu-Osher representation 

(7.3) 



f {i) dv = J2 (an, J 0f ij) dv + J <fyu ■ V x f ij) dv 

This moment equation is the same as the equation on p in the Euler system, and the classical 
proof for p being positive for the Euler equation can just be adopted [TT]. To check the positivity 
of T, one just need to make use of the last line of the moment equation, i.e. 



j 

i-l 



a 



f^dv + Pijh / 



i-l 



V x (fU) - dv 



(7.4a) 
(7.4b) 



(7.4a) is exactly what one could get when computing for E in the Euler system: the form of M 

2 

closes it up. So the classical method to prove that E > in Runge-Kutta scheme could be 
used, and the only thing new is from (7.4b). However, as proved in the section about AP, the 
difference between / and M is at most of e, thus (7.4b) is of order O(he). □ 
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7.2 |P(/) - P(g)\ <\f-g\ in d 2 norm 

We adopt the results from [28] . They denote P2 the collection of distributions F such that 



A metric d 2 on P2 is defined by 



/ \v\ 2 dF(v) < 00 



d 2 (F, G) = sup € 



(7.5) 



where / is the Fourier transform of F 

/(£) = J e-* v dF(v) 
One can transform the Boltzmann equation into its Fourier space and obtains [23} 12] 



/K + )/ ( D - / ( {)/(0) 



dn 



(7.6) 



where £ 



± _ g±lglg 



Theorem 5. d 2 (Pf,P g ) < d 2 {f,g) for Maxwell molecules with cut-off collision kernel. 
Proof. For Maxwell molecule with cut-off collision kernel J B = S. Thus 



sup|(5 I = sup 



Bf^dfldv* 



sup|p5| < 00. 



Considering P = Q + fj,f = Q + + (/x — Q ) /, it is enough to prove d 2 (Q^ , Q+) < Cd 2 (f,g) for 
C big enough. Given 



Q 



f 



B 



£ ■ n 

Iff 



dn, 



one has 



Qf-Qg 



B 



s 2 



£ ■ n 



f(e)f(n-m + )m- 



From |28J, one gets 



Thus, one has: 



\Z\ 2 

Q~f -Qj 



< 



sup 



dn 



f-g 



< S sup 



f-g 



l£l s 



Sd 2 (f,g) 



□ 
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